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Abstract A non-linear material model is proposed to 
describe the timber behavior. Anisotropic elasto-plastic 
constitutive law with hardening according to Hill yield cri¬ 
terion was used for the compressive behavior. Brittle 
behavior in tension and shear was modeled using the mod¬ 
ified Hill failure criterion. The material model was imple¬ 
mented in a finite element code to simulate embedding 
strength of glued laminated timber for dowel-type fasteners. 
Reasonable agreement is found between numerical simula¬ 
tions and experimental measurements. 

Keywords Dowel-type fasteners • Embedding strength • 
Finite element • Non-linearity 

Introduction 

Glued laminated timber, also called Glulam, is the most 
popular structural timber product in heavy wood con¬ 
structions. In fact, glulam is the only engineered wood 
product that can be produced in curved shapes, offering 
unlimited design flexibility. It is available in a range of 
appearance characteristics to meet end-use requirements. 
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In the timber structures, the connections are usually made 
with dowel-type fasteners and steel plates. 

According to the current European code for timber 
structures [1], the design of dowel-type timber connections 
is based on the European Yield Model proposed by 
Johansen [2]. Embedding strength (or dowel-bearing 
strength) is one of the key material property parameters for 
timber joint design based on this theory. Two embedding 
test methods are standardized such as in the EN383 stan¬ 
dard [3] and ASTM D 5764-97a [4]. Figure 1 shows the 
embedding test arrangements applied in the European and 
American methods. The European test method is based on 
the principle of reproducing the behavior of doweled 
connection but with short fastener to limit its bending 
deformations. Both tensile and compressive embedment 
testing modes are considered to take account of different 
stress distributions that prevail in real joints under lateral 
load. The clearance left between the fastener and the side- 
plates of the apparatuses minimizes the friction with cru¬ 
shed wood, but at the same time introduces the risks of 
increasing fastener deformations in bending during testing 
[5]. The American test method applies uniform contact 
along the fastener length during loading. 

In addition to the test method, the definition of embed¬ 
ding strength is different in European and American stan¬ 
dard. The embedding strength in the European method is 
calculated from the maximum load attained or the load 
corresponding to 5 mm of embedding displacement when 
the maximum load is not attained. The American method 
uses a 5 % diameter offset load to calculate the embedding 
strength. A comparison of the European and American 
procedures has been presented by Hilson and Pope [6]. 

Embedding strength is non-unique for any given speci¬ 
men of timber but is dependent upon several factors, such 
as size, cross-sectional shape and surface condition of 
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Fig. 1 Schematic illustrations 
of embedding test from EC (left) 
and ASTM (right) 
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fastener, orientation of fastener with respect to timber 
grain, degree of initial contact between fastener and wood 
[5]. To study all possible variations of those factors, the 
experimental way would be very time consuming. There is, 
therefore, a need to develop a numerical approach to 
investigate embedding strength more effectively [7]. 

In the last decade, some 2D finite element (FE) models 
have been proposed [8-10]. These models can be efficient 
to simulate the embedding behavior with thin wood 
members (such as the American test). However, they are 
limited with long fasteners such as dowel-type connections 
with thick timber members and slender dowels, which are 
three-dimensional problems (non-uniform stress distribu¬ 
tions across the thickness of members). Therefore, it is 
necessary to develop a 3D FE model to simulate the 
embedding behavior of timber with dowel-type fastener, in 
particular, for multi-fastener joints with thick timber 
members and slender dowels. The 3D FE models have the 
ability to predict the load carrying capacity of dowel timber 
connections as well as the 3D stress field. 

The wood is anisotropic material, the compressive and 
tensile strengths are different even in the same direction 
such as the parallel or perpendicular to grain directions. 
The wood presents ductile compressive behavior with 
softening in the direction parallel to grain and with hard¬ 
ening in the direction perpendicular to grain [11], and 
brittle failure in tension and shear. Thus, an appropriate 
wood constitutive law is important for the simulation of the 
non-linear behavior of embedding characteristic. Currently, 
wood is popularly modeled as elasto-plastic orthotropic in 
compression according to the Hill yield criterion. For 
simplification, the elastic perfectly plastic behavior has 
been adopted in some models [12-15]. However, the 
accurate prediction requires considering material hardening 
or softening [16, 17]. 

The progressive failure of timber occurs in the loading 
process resulting from brittle behavior in tension and shear. 
Few models have taken into account the contribution of 
tension or shear in the wood failure [10, 18-20]. It has long 


been suspected that stresses in principal directions influ¬ 
ence each other and therefore affect material failure. 
Especially for fibrous materials, stress interactions may 
change failure mode and failure stress. Multiaxial stresses 
must be considered in two- or three-dimensional structural 
elements such as plates, shells or elements under combined 
loading [21]. Some numerical models including the Tsai- 
Wu failure criterion, which is based on interactive failure 
theory, have been developed [13, 22]. A concern exists 
about using the Tsai-Wu criterion due to the need to 
evaluate the stress interaction coefficients that are sensitive 
to the test methods used to evaluate biaxial strength. No 
recommendations exist for selecting interaction coeffi¬ 
cients for 3D wood failure [23]. Although several failure 
criteria exist for wood and orthotropic materials, most are 
difficult to apply to 3D stress fields. 

It should be noted that few 3D FE models representing 
the behavior of timber connections are available in the 
literature. This paper represents a contribution in this field. 
A suitable material model for wood will allow engineer to 
predict embedding behavior in numerous applications 
without the reliance on extensive empirical tests. In this 
study, timber in compression is modeled as elasto-plastic 
material with hardening according to Hill yield criterion, 
and brittle behavior in tension and shear is modeled using 
the modified Hill failure criterion. Good correlation has 
been obtained between the experimental results and the FE 
simulations. 

Constitutive modeling 

Elasto-plastic behavior of wood in compression 

The timber properties are orthotropic with identical prop¬ 
erties in radial and tangential directions. This combined 
direction is referred to as perpendicular to grain (_L) while 
the longitudinal direction is referred to as parallel to grain 
(II). Wood has a cellular structure that is made up of 
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lumens and cell walls. Due to this structure, wood cells 
buckle and collapse under compressive loads. In com¬ 
pression parallel to grain, clear wood shows plastic soft- 


version of the von Mises yield criterion to take into account 
the anisotropy of the materials. Its stress potential can be 
expressed as follows: 
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ening, in compression perpendicular to grain, it shows 
ductile behavior with extended yielding plateau followed 
by significant hardening after the cell walls had collapsed 
into the lumens and the wood densified [11]. 

However, for the glued laminated timber, this plastic 
softening phenomenon in the longitudinal direction is not 
significant, it can be approximately regarded as elastic 
perfectly plastic [24]. In the direction perpendicular to 
grain, the behavior of glulam in compression can be rep¬ 
resented by a bi-linear elasto-plastic law including hard¬ 
ening [24]. The simplified bi-linear stress-strain curve is 
shown in Fig. 2, where the hardening slope (tan 6) is equal 
to 10 % of E 90 . 

The theory of plasticity is a procedure by a set of con¬ 
stitutive equations for a multi-axial stress state. This is 
accomplished on three basic properties: a yield criterion, a 
flow rule, and a hardening rule. Hill’s yield criterion (stress 
potentials) [25] has been adopted to simulate anisotropic 
plastic behavior of the glulam. This criterion is a generalized 



Fig. 2 Simplified bi-linear stress (cr)-strain(e) curve in compression 
perpendicular to grain 


The hypothesis of associated plasticity is considered. 
Thus, the Hill yield criterion is used as a plastic flow law. 
The relationship between the plastic strain increment and 
the stress increment is given by: 


d£ p = dz — (2) 

p da J 

where dA is the plastic multiplier. In this study, the 
isotropic hardening rule was adopted as in [17, 26], which 
assumes that the center of the yield surface remains 
stationary in the stress space, but that the size (radius) of 
the yield surface expands, due to hardening. 





H is the plastic modulus of the effective stress-effective 
strain diagram, which can be derived with the following set 
of equations [10]: 


= (4) 

E„ = ££t 
p E-Ej 

k 0 and k are the initial and updated effective yield 
stresses, respectively. E p is the plastic modulus, and E T is 
the tangent modulus of the stress-strain curve. Depending 
on the chosen value of //, perfect plasticity or hardening 
can be modeled. 



Brittle failure behavior in tension and shear 

In shear and tension, timber essentially exhibits behavior 
with brittle character which is represented by the modified 
Hill failure criterion. The Hill failure criterion was refor¬ 
mulated in such a way that tension and shear stresses which 
lead to brittle failure modes were considered. It can be 
expressed by Eq. (6). 
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Where and i/ 7 - are the stresses in timber, 

/t o and / 90 are the tensile strength parallel and per¬ 
pendicular to grain, 

f c 0 and f c 90 are the compressive strength parallel and 
perpendicular to grain, 

/ v is the shear strength of timber. 

Numerically, in this criterion, the progressive failure is 
simulated through a reduction of the material modulus to 
represent the damage evolution in timber. Thus, when the 
criterion is reached at an integration point, the elastic 
modulus in the direction parallel to grain E 0 and the 
modulus in the direction perpendicular to grain E 90 , are set 
equal to 10 % of their initial values. Hence, a kind of 
damage evolution in timber is considered. 


Table 1 Dimensions of specimens 


«(°) 

d (mm) 

L (mm) 

h (mm) 

t (mm) 

0 

16 

96 

112 

40 

45 

16 

160 

88 

40 

90 

16 

224 

64 

40 



Numerical applications 


Fig. 4 Finite element mesh for the embedding test 


Embedding test 

To validate the developed FE model, embedding tests were 
carried out with various orientations of fastener with respect 
to timber grain according to ASTM standard. The config¬ 
uration and geometry of the specimens are shown in Fig. 3 
and Table 1. 

Model description 


Table 2 Timber properties used in the model 


Eq (MPa) 

Ego (MPa) 

Gmean (MPa) 

V X R — V LT 

Vrl 

11600 

390 

750 

0.41 

0.02 


8 -noded hexahedral elements, is shown in Fig. 4. The 
model includes three parts: timber member, stiff steel 
dowel and plate at support. 


3D FE model is developed using the MARC.MSC software 
package [27]. The geometry of the FE model, based on 



Fig. 3 Configuration of embedding test 


Material properties 

The glulam used in the tests corresponds to the resistance 
class GF28h in accordance with standard “EN 1194”. Its 

o 

measured mean density p mean is equal to 453 kg/m and the 
moisture content /z mean is equal to 11 %. The values of the 
mechanical elastic parameters used are given in Table 2. 

In this study, timber tensile strength (parallel and per¬ 
pendicular to grain) and shear strength values are consid¬ 
ered: /t o = 19.8, / 9 Q = 0.48 and / v = 3.24 MPa. These 
mean values are obtained from the characteristic values 
using the coefficient 1.2. Considering a normal distribution, 
this equivalence corresponds approximately to a coefficient 
of variation equal to 10 %. Mean values of compressive 
strength are determined by test according to European 
standard: / c 0 = 39 and / c?90 = 5.74 MPa. As the defor¬ 
mation on the dowel is within the elastic range and also 
very small, therefore, the dowel is considered as an elastic 
material. 
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Load, boundary and contact conditions 

The loads are introduced using controlled displacements 
applied on the nodes of the steel dowel with a management 
of the contact condition between the steel dowel and the 
timber member. Contact was modeled using the direct 
constraint method in MSC.Marc. The method requires the 
definition of the “contact body” that potentially may come 
in contact with the other body. Contact bodies can simply 
be the physical bodies themselves (e.g. timber, dowel and 
steel plate). In this method, when a node of a body contacts 
another body, a multipoint constraint is imposed. This 
constraint allows the contacting node to slide on the con¬ 
tacted segment. In this way, a contacting node is forced to 
be on the contacted segment. During the iteration proce¬ 
dure, a node can slide from one segment to another, 
changing the retained nodes associated with the constraint. 
A node is considered sliding off a contacted segment if it 
passes the end of the segment over a distance more than the 
contact tolerance. 

Contact may be developed with friction based on the 
Coulomb criterion. The method allows no movement until 
the friction force is reached. After that, the movement is 
initiated and the friction forces remain constant. The fric¬ 
tion coefficient between the timber member and stiff steel 
plate was set equal to 0.2 [28]. Considering surface con¬ 
dition of dowel, the friction coefficient between the dowel 
and timber was set equal to 0.3, which can be considered as 
the limit value between the smooth and rough surface 
dowel [29]. 

Model validation 

The 3D FE models are validated on the basis of a com¬ 
parison between the embedding strength-displacement 
curves given by FE models and those from tests. Three 
series of experiments were conducted, and four specimens 
were tested in each series. The specimens were labeled as 



Fig. 5 Comparison between numerical and experimental embedding 
strength-displacement curves under loading parallel to grain 


fh-x-a,b,c,d: “x” represented the angle between the loading 
direction and the direction of grain; “a,b,c,d” represented 
the specimen identification. The predicted embedding 
strength-displacement curve under loading parallel to grain 
is presented in Fig. 5. Though the predicted embedding 
strength is lower than the experimental ones, the model 
shows a similar tendency with the experimental curves. 
The underestimate of the embedding strength can be 
explained by the fact that the FE model does not take into 
account the material densification effect underneath the 
dowel. In fact, due to densification, wood underneath the 
dowel should have higher density and modulus. 

Figures 6 and 7 show the reasonable prediction of the 
embedding strength-displacement curves under loading at 
45° from the longitudinal direction and under loading 
perpendicular to grain. The differences between numerical 
and experimental curves, on the hardening slopes in the 
initial plasticity range, can be explained by the simplifi¬ 
cation of bi-linear compression behavior with hardening in 
the direction perpendicular to grain. In fact, after a linear 
elastic part, the crushing strength of wood fibers is reached. 
Then, a certain kind of plateau occurs and, at higher strains, 
the curve increases as the local density of wood increases 
(Fig. 2). 



Fig. 6 Comparison between numerical and experimental embedding 
strength-displacement curves under the loading at 45° from the 
longitudinal direction 



Fig. 7 Comparison between numerical and experimental embedding 
strength-displacement curves under loading perpendicular to grain 
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a 90 (MPa) 

0.480 

0.432 

0.384 

0.336 

0.288 

0.240 

0.192 

0.144 

0.096 

0.048 

0.000 




Fig. 8 Perpendicular to grain tensile stress distribution at 5 % 
diameter offset load parallel to grain 


I (MPa) 

3.240 

2.592 

1.944 

1.296 

0.648 

0.000 

- 0.648 

- 1.296 

- 1.944 

- 2.592 

- 3.240 




Load direction 0° 



Load direction 45° 



Load direction 90° 


Fig. 9 Shear stress distribution at 5 % diameter offset load parallel to 
grain 

Stress analysis 

For tests with specimens loaded in the direction parallel to 
grain, the failure occurred prematurely by splitting parallel 
to grain. The brittle failure is related to tensile stresses 
perpendicular to grain and shear stresses. Figures 8 and 9 
present the tensile and shear stress distributions in timber 
at 5 % diameter offset load parallel to grain. The crack 


Fig. 10 Contour plot of the zones where the failure index is reached 

parallel to grain appeared and propagated under the dowel 
during the test, which is in accordance with the observed 
stress contours. 

The failure index of the modified Hill criterion based on 
the stresses’ interaction provided by the numerical model 
is used to evaluate the potential brittle failure position. 
Figure 10 shows the zone where the failure criterion is 
reached at 5 % diameter offset load. Under loading parallel 
to grain, the zone with high values of failure index is 


<£) Springer 






























































































































































































































































































































































































































































































J Wood Sci (2013) 59:17-23 


23 


located under dowel hole, which is consistent with the 
location of initial crack. For loading at 45° from the lon¬ 
gitudinal direction, the shear failure is close to the dowel 
hole and along the direction parallel to grain, which is 
consistent with the numerical observation. Under loading 
perpendicular to grain, the main failure is the embedding 
failure, due to the crushing of wood fibers, the splitting 
parallel to grain occurred on both sides of the dowel hole. 
The brittle failure can be predicted by this model. 

Conclusion 

A 3D FE model for glued laminated timber behavior is 
presented and successfully applied to predict the dowel 
type embedding strength. The proposed model describes 
the material non-linearity of timber in terms of plastic 
behavior in compression and brittle failure in tension and 
shear. Particularly for compression perpendicular to grain, 
the model can be used to simulate the strain hardening of 
timber. The developed model is an improvement on pre¬ 
vious wood constitutive models. It can be helpful to design 
and analyze various timber connections. 
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